library(reticulate)
use_condaenv("glmhmm")
source_python("python/glm_hmm.py")
import numpy as np
import numpy.random as npr
import matplotlib.pyplot as plt
import ssm
quit
library(tidyverse)
library(broom)
library(dplyr)
library(R.matlab)
library(cowplot)
list.files("R", full.names = TRUE) |>
walk(source)
load("RData/glm_hmm1_seeds.RData")
load("RData/glm_hmm2_seeds.RData")
load("RData/glm_hmm3_seeds.RData")
Renames the states according to their values (so similar states
between subjects can be identified)
# glm_hmm1 <- glm_hmm1_seeds %>%
# mutate(fit = list(glm_hmm_extract_best(fit))) %>%
# order_states_glm_hmm()
#
# glm_hmm1
#
# save(glm_hmm1, file = "RData/glm_hmm1.RData")
# glm_hmm2 <- glm_hmm2_seeds %>%
# mutate(fit = list(glm_hmm_extract_best(fit))) %>%
# order_states_glm_hmm()
#
# glm_hmm2
#
# save(glm_hmm2, file = "RData/glm_hmm2.RData")
# glm_hmm3 <- glm_hmm3_seeds %>%
# mutate(fit = list(glm_hmm_extract_best(fit))) %>%
# order_states_glm_hmm()
#
# glm_hmm3
#
# save(glm_hmm3, file = "RData/glm_hmm3.RData")
load("RData/glm_hmm1.RData")
load("RData/glm_hmm2.RData")
load("RData/glm_hmm3.RData")
———————————————————————
MODEL COMPARISON
liks de cada modelo
DL:
liks_glm_hmm <- liks_glm_hmm1 |>
bind_rows(liks_glm_hmm2) |>
bind_rows(liks_glm_hmm3) |>
mutate(aic = -2 * log_lik + 2 *n_par)
save(liks_glm_hmm, file = "RData/liks_glm_hmm.RData")
liks_glm_hmm_best <- liks_glm_hmm |>
group_by(subject) |>
filter(aic == min(aic))
liks_glm_hmm_best |>
ungroup() |>
count(model)
save(liks_glm_hmm_best, file = "RData/liks_glm_hmm_best.RData")
LS0tCnRpdGxlOiAiaGxtaG1tIG1ldGFjb2duaXRpb24gc2VlZCBzZWxlY3Rpb24gYW5kIG1vZGVsIGNvbXBhcmlzb24iCm91dHB1dDogaHRtbF9ub3RlYm9vawotLS0KCgpgYGB7cn0KbGlicmFyeShyZXRpY3VsYXRlKQp1c2VfY29uZGFlbnYoImdsbWhtbSIpCgpzb3VyY2VfcHl0aG9uKCJweXRob24vZ2xtX2htbS5weSIpCmBgYAoKYGBge3B5dGhvbn0KaW1wb3J0IG51bXB5IGFzIG5wCmltcG9ydCBudW1weS5yYW5kb20gYXMgbnByCmltcG9ydCBtYXRwbG90bGliLnB5cGxvdCBhcyBwbHQKaW1wb3J0IHNzbSAKYGBgCgpgYGB7cn0KbGlicmFyeSh0aWR5dmVyc2UpCmxpYnJhcnkoYnJvb20pCmxpYnJhcnkoZHBseXIpCmxpYnJhcnkoUi5tYXRsYWIpCmxpYnJhcnkoY293cGxvdCkKCmxpc3QuZmlsZXMoIlIiLCBmdWxsLm5hbWVzID0gVFJVRSkgfD4gCiAgd2Fsayhzb3VyY2UpCmBgYAoKCmBgYHtyfQpsb2FkKCJSRGF0YS9nbG1faG1tMV9zZWVkcy5SRGF0YSIpCmxvYWQoIlJEYXRhL2dsbV9obW0yX3NlZWRzLlJEYXRhIikKbG9hZCgiUkRhdGEvZ2xtX2htbTNfc2VlZHMuUkRhdGEiKQpgYGAKCgpSZW5hbWVzIHRoZSBzdGF0ZXMgYWNjb3JkaW5nIHRvIHRoZWlyIHZhbHVlcyAoc28gc2ltaWxhciBzdGF0ZXMgYmV0d2VlbiBzdWJqZWN0cyBjYW4gYmUgaWRlbnRpZmllZCkKYGBge3J9CiMgZ2xtX2htbTEgPC0gZ2xtX2htbTFfc2VlZHMgJT4lIAojICAgICBtdXRhdGUoZml0ID0gbGlzdChnbG1faG1tX2V4dHJhY3RfYmVzdChmaXQpKSkgJT4lCiMgICAgIG9yZGVyX3N0YXRlc19nbG1faG1tKCkKIyAKIyBnbG1faG1tMQojIAojIHNhdmUoZ2xtX2htbTEsIGZpbGUgPSAiUkRhdGEvZ2xtX2htbTEuUkRhdGEiKQpgYGAKCmBgYHtyfQojIGdsbV9obW0yIDwtIGdsbV9obW0yX3NlZWRzICU+JSAKIyAgIG11dGF0ZShmaXQgPSBsaXN0KGdsbV9obW1fZXh0cmFjdF9iZXN0KGZpdCkpKSAlPiUKIyAgIG9yZGVyX3N0YXRlc19nbG1faG1tKCkKIyAKIyBnbG1faG1tMgojIAojIHNhdmUoZ2xtX2htbTIsIGZpbGUgPSAiUkRhdGEvZ2xtX2htbTIuUkRhdGEiKQpgYGAKCmBgYHtyfQojIGdsbV9obW0zIDwtIGdsbV9obW0zX3NlZWRzICU+JQojICAgbXV0YXRlKGZpdCA9IGxpc3QoZ2xtX2htbV9leHRyYWN0X2Jlc3QoZml0KSkpICU+JQojICAgb3JkZXJfc3RhdGVzX2dsbV9obW0oKQojIAojIGdsbV9obW0zCiMgCiMgc2F2ZShnbG1faG1tMywgZmlsZSA9ICJSRGF0YS9nbG1faG1tMy5SRGF0YSIpCmBgYAoKYGBge3J9CmxvYWQoIlJEYXRhL2dsbV9obW0xLlJEYXRhIikKbG9hZCgiUkRhdGEvZ2xtX2htbTIuUkRhdGEiKQpsb2FkKCJSRGF0YS9nbG1faG1tMy5SRGF0YSIpCmBgYAoKCi0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLQotLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0KTU9ERUwgQ09NUEFSSVNPTiAKLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tCi0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLQoKCmxpa3MgZGUgY2FkYSBtb2RlbG8KYGBge3J9Cmxpa3NfZ2xtX2htbTEgPC0gZ2xtX2htbTFfc2VlZHMgJT4lIAogIG11dGF0ZShmaXQgPSBsaXN0KGdsbV9obW1fZXh0cmFjdF9iZXN0KGZpdCkpKSAlPiUgCiAgdW5uZXN0X3dpZGVyKGZpdCkgJT4lIAogIG11dGF0ZShuID0gbWFwKGRhdGEsIGNvdW50KSkgJT4lIAogIHVubmVzdChuKSAlPiUgCiAgc2VsZWN0KHN1YmplY3QsIGxvZ19saWssIG5fcGFyLCBuKSAlPiUgCiAgbXV0YXRlKG1vZGVsID0gImludC1zbG8gMSIpCgpsaWtzX2dsbV9obW0yIDwtIGdsbV9obW0yX3NlZWRzICU+JSAKICBtdXRhdGUoZml0ID0gbGlzdChnbG1faG1tX2V4dHJhY3RfYmVzdChmaXQpKSkgJT4lIAogIHVubmVzdF93aWRlcihmaXQpICU+JSAKICBtdXRhdGUobiA9IG1hcChkYXRhLCBjb3VudCkpICU+JSAKICB1bm5lc3QobikgJT4lIAogIHNlbGVjdChzdWJqZWN0LCBsb2dfbGlrLCBuX3BhciwgbikgJT4lIAogIG11dGF0ZShtb2RlbCA9ICJpbnQtc2xvIDIiKQoKbGlrc19nbG1faG1tMyA8LSBnbG1faG1tM19zZWVkcyAlPiUKICBtdXRhdGUoZml0ID0gbGlzdChnbG1faG1tX2V4dHJhY3RfYmVzdChmaXQpKSkgJT4lCiAgdW5uZXN0X3dpZGVyKGZpdCkgJT4lCiAgbXV0YXRlKG4gPSBtYXAoZGF0YSwgY291bnQpKSAlPiUKICB1bm5lc3QobikgJT4lCiAgc2VsZWN0KHN1YmplY3QsIGxvZ19saWssIG5fcGFyLCBuKSAlPiUKICBtdXRhdGUobW9kZWwgPSAiaW50LXNsbyAzIikKYGBgCgpETDoKYGBge3J9Cmxpa3NfZ2xtX2htbSA8LSBsaWtzX2dsbV9obW0xIHw+IAogIGJpbmRfcm93cyhsaWtzX2dsbV9obW0yKSB8PiAKICBiaW5kX3Jvd3MobGlrc19nbG1faG1tMykgfD4gCiAgbXV0YXRlKGFpYyA9IC0yICogbG9nX2xpayArIDIgKm5fcGFyKQoKc2F2ZShsaWtzX2dsbV9obW0sIGZpbGUgPSAiUkRhdGEvbGlrc19nbG1faG1tLlJEYXRhIikKCmxpa3NfZ2xtX2htbV9iZXN0IDwtIGxpa3NfZ2xtX2htbSB8PiAKICBncm91cF9ieShzdWJqZWN0KSB8PiAKICBmaWx0ZXIoYWljID09IG1pbihhaWMpKSAKCmxpa3NfZ2xtX2htbV9iZXN0IHw+IAogIHVuZ3JvdXAoKSB8PiAKICBjb3VudChtb2RlbCkKCnNhdmUobGlrc19nbG1faG1tX2Jlc3QsIGZpbGUgPSAiUkRhdGEvbGlrc19nbG1faG1tX2Jlc3QuUkRhdGEiKQpgYGAKCg==